{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "#from numpy.linalg import cholesky\n",
    "import matplotlib.pyplot as plt\n",
    "from functools import reduce"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "sampleNo = 100\n",
    "mu = np.array([[0, 2]])\n",
    "mu1 = np.array([[1.5, 0]])\n",
    "Sigma = np.array([[2, 1], [1, 2]])\n",
    "#R = cholesky(Sigma).T\n",
    "#va,vc = np.linalg.eig(Sigma); R2 = (np.diag(va)**0.5)@vc.T\n",
    "\n",
    "s1 = np.random.multivariate_normal(mu[0],Sigma,sampleNo) \n",
    "s2 = np.random.multivariate_normal(mu1[0],Sigma,sampleNo) \n",
    "\n",
    "x = range(-100,100,1)\n",
    "C = np.linalg.inv(Sigma)\n",
    "w = 2*np.dot(C, mu1.T-mu.T)\n",
    "b = reduce(np.dot, [mu,C,mu.T])- reduce(np.dot, [mu1,C,mu1.T])\n",
    "w0 = (-w[0]/w[1])[0]\n",
    "y= np.dot(w0, x)+(-b/w[1])\n",
    "\n",
    "plt.plot(*s1.T,'.',label = 's1')\n",
    "plt.plot(*s2.T,'.',label = 's2')\n",
    "plt.axis('scaled')\n",
    "plt.plot(mu[0],mu1[0],'-')\n",
    "plt.plot(x,y[0],'-')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
